MSK SPECTRUM freeze major subset deep dives

library(magrittr)
library(tidyverse)
library(Seurat)
library(readxl)
library(cowplot)
library(colorblindr)
library(viridis)
library(progeny)
library(destiny)

coi <- params$cell_type_super
cell_sort <- params$cell_sort
cell_type_major <- params$cell_type_major
louvain_resolution <- params$louvain_resolution
louvain_cluster <- params$louvain_cluster

1 Cluster markers

1.1 Major T.super markers for cell assign

### load all data ---------------------------------
source("_src/global_vars.R")

# seu_obj <- read_rds(paste0("/work/shah/isabl_data_lake/analyses/16/52/1652/celltypes/", coi, "_processed.rds"))
seu_obj <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/outs_pre/", coi, "_seurat_", louvain_resolution, ".rds"))
# seu_obj <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_highqc.rds"))

myfeatures <- c("umapharmony_1", "umapharmony_2", "sample", louvain_cluster, "doublet", "nCount_RNA", "nFeature_RNA", "percent.mt", "doublet_score")

plot_data_wrapper <- function(cluster_res) {
  cluster_res <- enquo(cluster_res)
  as_tibble(FetchData(seu_obj, myfeatures)) %>% 
    left_join(meta_tbl, by = "sample") %>% 
    rename(cluster = !!cluster_res) %>% 
    mutate(cluster = as.character(cluster),
           tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite)))
}

plot_data <- plot_data_wrapper(louvain_cluster)

1.2 Subtype currated markers

helper_f <- function(x) ifelse(is.na(x), "", x)

markers_v6_super[[coi]] %>% 
  group_by(subtype) %>% 
  mutate(rank = row_number(gene)) %>% 
  spread(subtype, gene) %>% 
  mutate_all(.funs = helper_f) %>% 
  formattable::formattable()

rank	CD4.T.CXCL13	CD4.T.naive	CD4.T.reg	CD4.Th17	CD8.T.CXCL13	CD8.T.cytotoxic	CD8.T.ISG	Cycling.T.NK	dissociated	NK.CD56	NK.cytotoxic
1	CD4	CCR7	CD4	CCR6	CCL3	CCL4	IFI6	ASPM	BTG1	AREG	FCGR3A
2	CD40LG	CD4	FOXP3	IL4I1	CCL4L2	CCL5	IFIT1	CENPF	DNAJB1	FCER1G	FGFBP2
3	CTLA4	CD40LG	IL2RA	IL7R	CD8A	CD8A	IFIT2	HIST1H4C	DUSP1	GNLY	KLRD1
4	CXCL13	IL7R	TNFRSF4	KLRB1	CRTAM	CD8B	IFIT3	HMGB2	EGR1	KLRC1	KLRF1
5	FKBP5	KLF2	TRAC	LST1	CXCL13	CRTAM	ISG15	MKI67	FOS	KRT81	PRF1
6	IL6ST	LTB		LTB	FABP5	CST7	MX1	STMN1	FOSB	TRDC	SPON2
7	ITM2A	TCF7		TNFSF13B	GZMB	DTHD1	MX2	TOP2A	HSPA1A	TYROBP
8	MAF	TPT1			HAVCR2	GZMA	RSAD2	TUBA1B	HSPA1B	XCL1
9	NMB				IFNG	GZMH		TUBB	HSPA6	XCL2
10	NR3C1				LAG3	GZMK		TYMS	JUN
11	PDCD1				MIR155HG	GZMM			JUNB
12	TNFRSF4				PHLDA1	HLA-DPB1			KLF2
13	TOX2				PTMS	ITM2C			MT1E
14	TSHZ2				RBPJ	KLRG1			MT1X
15					TNFRSF9	TRGC2

1.3 Subtype cluster markers

## define patient specific clusters
patient_specific_clusters <- plot_data %>% 
  group_by(patient_id_short, cluster) %>% 
  tally %>% 
  group_by(cluster) %>% 
  mutate(nrel = n / sum(n),
         ntotal = sum(n)) %>% 
  arrange(desc(nrel)) %>% 
  distinct(cluster, .keep_all = T) %>% 
  filter(nrel > 0.5) %>% 
  ungroup %>% 
  mutate(cluster_label_ps = make.names(paste0("OV.", patient_id_short, ".specific"), unique = T),
         cluster = as.numeric(cluster)) %>% 
  distinct(cluster, cluster_label_ps)



## Hypergeometric test --------------------------------------

# marker_tbl <- read_tsv(paste0("/work/shah/isabl_data_lake/analyses/16/52/1652/celltypes/", coi, "_markers.tsv")) %>% 
#   filter(resolution == louvain_resolution)
marker_tbl <- read_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/outs_pre/", coi, "_markers_", louvain_resolution, ".tsv"))
# marker_tbl <- read_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_highqc_markers_02.tsv"))

test_set <- marker_tbl %>% 
  group_by(cluster) %>% 
  filter(!str_detect(gene, "^RPS|^RPL")) %>% 
  slice(1:50) %>% 
  mutate(k = length(cluster)) %>% 
  ungroup %>%
  select(cluster, gene, k) %>% 
  mutate(join_helper = 1) %>% 
  group_by(cluster, join_helper, k) %>%
  nest(test_set = gene)

markers_doub_tbl <- markers_v6 %>% 
  enframe("subtype", "gene") %>% 
  filter(!(subtype %in% unique(c(coi, cell_type_major)))) %>% 
  unnest(gene) %>% 
  group_by(gene) %>% 
  filter(length(gene) == 1) %>% 
  mutate(subtype = paste0("doublet.", subtype)) %>% 
  bind_rows(tibble(subtype = "Mito.high", gene = grep("^MT-", rownames(seu_obj), value = T)))

ref_set <- markers_v6_super[[coi]] %>% 
  bind_rows(markers_doub_tbl) %>% 
  group_by(subtype) %>% 
  mutate(m = length(gene),
         n = length(rownames(seu_obj))-m,
         join_helper = 1) %>% 
  group_by(subtype, m, n, join_helper) %>%
  nest(ref_set = gene)

hyper_tbl <- test_set %>% 
  left_join(ref_set, by = "join_helper") %>% 
  group_by(cluster, subtype, m, n, k) %>%
  do(q = length(intersect(unlist(.$ref_set), unlist(.$test_set)))) %>%
  mutate(pval = 1-phyper(q = q, m = m, n = n, k = k)) %>%
  ungroup %>%
  mutate(qval = p.adjust(pval, "BH"),
         sig = qval < 0.01)

# hyper_tbl %>% 
#   group_by(subtype) %>% 
#   filter(any(qval < 0.01)) %>%
#   ggplot(aes(subtype, -log10(qval), fill = sig)) +
#   geom_bar(stat = "identity") +
#   facet_wrap(~cluster) +
#   coord_flip()
  
low_rank <- str_detect(unique(hyper_tbl$subtype), "doublet|dissociated")
subtype_lvl <- c(sort(unique(hyper_tbl$subtype)[!low_rank]), sort(unique(hyper_tbl$subtype)[low_rank]))
  
cluster_label_tbl <- hyper_tbl %>% 
  mutate(subtype = ordered(subtype, levels = subtype_lvl)) %>% 
  arrange(qval, subtype) %>%
  group_by(cluster) %>% 
  slice(1) %>% 
  mutate(subtype = ifelse(sig, as.character(subtype), paste0("unknown_", cluster))) %>% 
  select(cluster, cluster_label = subtype) %>% 
  ungroup %>% 
  mutate(cluster_label = make.unique(cluster_label, sep = "_")) %>% 
  left_join(patient_specific_clusters, by = "cluster") %>% 
  mutate(cluster_label = ifelse(is.na(cluster_label_ps), cluster_label, cluster_label_ps)) %>% 
  select(-cluster_label_ps)

seu_obj$cluster_label <- unname(deframe(cluster_label_tbl)[as.character(unlist(seu_obj[[paste0("RNA_snn_res.", louvain_resolution)]]))])
plot_data$cluster_label <- seu_obj$cluster_label

cluster_n_tbl <- seu_obj$cluster_label %>% 
  table() %>% 
  enframe("cluster_label", "cluster_n") %>% 
  mutate(cluster_nrel = cluster_n/sum(cluster_n))

marker_sheet <- marker_tbl %>% 
  left_join(cluster_label_tbl, by = "cluster") %>% 
  group_by(cluster_label) %>% 
  filter(!str_detect(gene, "^RPS|^RPL")) %>% 
  slice(1:50) %>% 
  mutate(rank = row_number(-avg_logFC)) %>% 
  select(cluster_label, gene, rank) %>% 
  ungroup %>% 
  mutate(cluster_label = ordered(cluster_label, levels = unique(c(names(clrs$cluster_label[[coi]]), sort(cluster_label))))) %>% 
  spread(cluster_label, gene) %>% 
  mutate_all(.funs = helper_f)

marker_tbl_annotated <- marker_tbl %>% 
  left_join(cluster_label_tbl, by = "cluster") %>% 
  left_join(cluster_n_tbl, by = "cluster_label") %>% 
  select(-cluster, -resolution) %>% 
  mutate(cluster_label = ordered(cluster_label, levels = unique(c(names(clrs$cluster_label[[coi]]), sort(cluster_label))))) %>% 
  arrange(cluster_label, -avg_logFC, p_val_adj)

write_tsv(marker_sheet, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_sheet.tsv"))

write_tsv(marker_sheet, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_sheet.tsv"))

write_tsv(marker_tbl_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated.tsv"))

write_tsv(marker_tbl_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_table_annotated.tsv"))

formattable::formattable(marker_sheet)

rank	CD4.T.naive	CD4.T.CXCL13	CD4.T.reg	CD4.Th17	CD8.T.cytotoxic	CD8.T.CXCL13	CD8.T.ISG	NK.CD56	NK.cytotoxic	Cycling.T.NK	Cycling.T.NK_1	dissociated	dissociated_1	doublet.Fibroblast	doublet.Monocyte	OV.008.specific	OV.008.specific.1	OV.085.specific
1	IL7R	CXCL13	TNFRSF4	KLRB1	GZMK	CXCL13	IFIT3	GNLY	FGFBP2	STMN1	CENPF	CCL4L2	HSPA1A	DCN	HLA-DRA	KLF2	SOX4	IGFBP5
2	CCR7	NMB	IL2RA	IL4I1	CD8A	GZMB	ISG15	TYROBP	FCGR3A	TUBA1B	ASPM	CCL4	HSPA1B	IGFBP5	CST3	CCR7	PTCRA	TAGLN
3	KLF2	NR3C1	FOXP3	IL7R	CD8B	CCL4L2	MX1	AREG	SPON2	MKI67	MKI67	IFNG	MT1X	RBP1	CXCL8	JUNB	MAL	ADIRF
4	EEF1B2	MAF	CTLA4	LTB	ITM2C	MIR155HG	IFIT1	KLRC1	PRF1	HIST1H4C	HMGB2	FOS	DNAJB1	C7	SPP1	FOS	MZB1	DCN
5	TPT1	FKBP5	LTB	LST1	GZMH	TNFRSF9	RSAD2	FCER1G	KLRF1	TUBB	TOP2A	FOSB	MT1E	MEG3	S100A9	DUSP1	DNTT	IGFBP7
6	EEF1A1	IL6ST	RTKN2	TNFSF13B	CCL5	HAVCR2	IFIT2	TRDC	GNLY	TOP2A	UBE2C	TNF	HSPA6	CALD1	S100A8	SELL	TFDP2	SPARCL1
7	MAL	ITM2A	BATF	CCR6	TRGC2	RBPJ	IFI6	XCL1	KLRD1	TYMS	CCNB1	JUN	HSP90AA1	IGFBP4	HLA-DQA1	IL7R	CD1E	CALD1
8	TCF7	TSHZ2	TNFRSF18	CTSH	KLRG1	LAG3	MX2	KLRD1	NKG7	CENPF	UBE2S	EGR1	HSPH1	RARRES2	CD74	AREG	STMN1	MGP
9	CD40LG	CTLA4	SAT1	AQP3	GZMA	IFNG	IFI44L	KRT81	CX3CR1	HMGB2	PTTG1	CCL3L1	HSPE1	NR2F2	LYZ	EEF1A1	ARPP21	C11orf96
10	SELL	CD40LG	TBC1D4	CCL20	CCL4	CCL3	ISG20	XCL2	GZMB	ASPM	TPX2	NFKBID	HSPD1	SELENOP	FTL	CD69	AC084033.3	MYL9
11	GPR183	PDCD1	TIGIT	NFKBIA	CRTAM	PTMS	HERC5	IGFBP2	PLAC8	NUSAP1	STMN1	CD69	HSPB1	MDK	APOE	BTG2	CDK6	TIMP3
12	LDHB	CD4	GADD45A	RORA	CST7	CD8A	SAMD9L	CLIC3	CLIC3	PCLAF	KPNA2	NR4A2	HSP90AB1	SOX4	HLA-DQB1	GPR183	MAP1A	MDK
13	NOSIP	LIMS1	TNFRSF1B	TNFRSF25	GZMM	CRTAM	OAS1	IL2RB	PLEK	HMGN2	CENPE	AC020916.1	DNAJA1	ADIRF	MARCKS	PIK3IP1	AC011893.1	TPM2
14	SNHG8	TNFRSF4	PMAIP1	CEBPD	HLA-DPB1	PHLDA1	TNFSF10	CEBPD	TYROBP	H2AFZ	CKS2	EGR2	JUN	MGP	AIF1	CD55	GLUL	NR2F2
15	PABPC1	RNF19A	UGP2	TNFAIP3	DTHD1	FABP5	STAT1	KRT86	PTGDS	PCNA	TUBA1B	DUSP1	CACYBP	EGR1	C1QB	DNAJB1	ADA	FBLN1
16	NOP53	CORO1B	IKZF2	NCR3	PPP1R14B	TIGIT	EIF2AK2	TXK	EFHD2	HIST1H1B	HMMR	TNFSF9	HSPA8	STAR	FTH1	FKBP5	AL357060.1	IGF1
17	LEF1	RBPJ	TNFRSF9	SLC4A10	CD3G	KRT86	OAS3	CTSW	FCER1G	DUT	DLGAP5	KLF6	FKBP4	SFRP4	BASP1	TSC22D3	GRASP	RAMP1
18	EIF3E	CPM	ICOS	TMIGD2	THEMIS	JAML	SAMD9	KLRB1	CST7	CLSPN	CCNB2	TNFAIP3	CHORDC1	CLU	C1QA	RACK1	CD1B	SELENOP
19	LTB	ZBED2	LINC01943	DPP4	DUSP2	CCL5	XAF1	MATK	GZMH	SMC4	TUBB4B	DUSP2	RGS2	ADAMTS1	C15orf48	LDHB	MIR181A1HG	RARRES2
20	RACK1	AC004585.1	IL32	TPT1	CD3D	LINC01871	GBP1	CCL3	ADGRG1	ATAD2	NUSAP1	BTG2	FOS	TIMP2	FN1	SARAF	CCDC26	LUM
21	NACA	TOX2	SOX4	MYBL1	LYAR	CXCR6	EPSTI1	CD7	CCL3	SMC2	BIRC5	IER2	ANXA1	TCEAL4	MNDA	CXCR4	JCHAIN	IGFBP6
22	UBA52	DUSP4	ARID5B	S100A4	TC2N	TNIP3	IFI44	NKG7	HOPX	TPX2	HMGN2	PPP1R15A	DNAJB4	C1R	G0S2	PLAC8	VIPR2	CARMN
23	TOMM7	AHI1	CD27	CD40LG	EOMES	PDCD1	PLSCR1	CD63	IGFBP7	TMPO	ARL6IP1	NR4A1	PPP1R15A	WFDC2	APOC1	EEF1B2	ID1	IGFBP4
24	SOCS3	ICA1	BIRC3	SPOCK2	CXCR6	HLA-DRB1	IFI35	HOPX	ZEB2	HELLS	CDC20	NFKBIA	ZFAND2A	FHL2	SOD2	TCF7	SOCS2	DST
25	JUNB	ARID5B	LAYN	LINC01871	HLA-DPA1	GAPDH	USP18	TMIGD2	PRSS23	UBE2C	CDKN3	JUND	FOSB	IFITM3	NPC2	TPT1	RCAN1	CAV1
26	SERINC5	CD84	CORO1B	ERN1	HLA-DRB1	FAM3C	OASL	CMC1	AKR1C3	NASP	SMC4	GADD45B	DNAJA4	PEG3	LST1	FOSB	CLDN5	ID4
27	TMEM123	CCDC50	TYMP	JAML	SLF1	CTLA4	MT2A	TNFRSF18	CD247	DEK	CKS1B	TSC22D3	DUSP1	C1S	GSN	NOSIP	CASC15	NUPR1
28	EEF2	IGFL2	DUSP4	FKBP11	APOBEC3G	SPRY1	HELZ2	KLRC2	MYBL1	RRM2	KIF20B	TAGAP	TSPYL2	LUM	MS4A6A	SC5D	PFKFB2	CSRP2
29	FXYD5	RGS1	CD4	IFNGR1	PPP2R5C	CCND2	TRIM22	SRGAP3	AREG	CKS1B	GTSE1	ZFP36	SERPINH1	SERPINF1	IL1B	NACA	GALNT2	CDKN1C
30	TSHZ2	BATF	ENTPD1	MGAT4A	KIAA1551	CD63	CMPK2	GSTP1	C1orf21	KNL1	TUBA1C	ZFP36L1	UBC	CEBPD	PSAP	AP3M2	HES4	PGR
31	TRABD2A	SRGN	CTSC	S100A6	F2R	CD8B	LY6E	LAT2	S1PR5	HMGB1	TUBB	RGCC	UBB	AKAP12	GRN	BTG1	MARCKSL1	COL6A1
32	ANK3	CH25H	MIR4435-2HG	ELK3	SLAMF7	GZMH	PARP14	GZMB	CTSW	MCM7	SGO2	IL7R	AHSA1	TIMP1	CD83	NOP53	TP53INP1	COL6A2
33	SARAF	SPOCK2	LINC02099	EEF1A1	CXCR4	ENTPD1	NT5C3A	LINC00996	ABHD17A	FABP5	H2AFZ	CRTAM	NEU1	CST3	CTSH	ZBTB16	APBA2	EMX2
34	AQP3	ZNRF1	MAGEH1	KIT	SH2D1A	SRGAP3	DDX58	PRF1	KLF2	UBE2S	JPT1	KDM6B	DNAJB6	NR2F1	GLUL	ERAP2	TSHR	PALLD
35	RIPOR2	CHN1	SPOCK2	LTC4S	FAM102A	VCAM1	IRF7	KLRF1	PTPN12	GAPDH	CEP55	NR4A3	GADD45B	DLK1	MEF2C	CCND3	NREP	PPP1R14A
36	AP3M2	TNFRSF25	PHACTR2	RUNX2	CD52	ID2	RNF213	NCAM1	TTC38	CXCL13	NUF2	ATF3	KLF6	SERPING1	CYBB	EEF1D	MME	RBP1
37	ZFAS1	CD200	CARD16	RORC	YBX3	HLA-DRA	DDX60L	CXXC5	KLRB1	RANBP1	KIF14	DDX3X	BTG2	BEX3	CTSB	PPP1R15A	AC002454.1	C7
38	LINC02273	METTL8	S100A4	ZBTB16	STK17A	ITGAE	DDX60	SH2D1B	CCL4	H2AFX	KNL1	DUSP6	JUNB	C11orf96	CSF3R	PLK3	CHI3L2	MFGE8
39	EIF4B	RILPL2	STAM	FAM241A	CCR5	DUSP4	SAT1	MCTP2	PTGDR	MCM3	CENPA	GZMK	TNF	FILIP1L	CD14	ZFP36L2	SMIM3	KANK2
40	TOB1	TNFRSF18	GLRX	PDE4D	GPR174	LYST	PPM1K	IFITM3	ITGB2	EZH2	CDK1	KLRG1	ERN1	RARRES1	C1QC	EEF2	SSBP2	NR2F1
41	SESN3	SLA	SPATS2L	IL23R	COTL1	TNFSF4	PNPT1	ZNF683	XBP1	TUBB4B	HMGB3	JUNB	JUND	COL1A2	SGK1	HNRNPA1	UHRF1	MFAP4
42	NSA2	SMCO4	AC005224.3	B3GALT2	CCL4L2	NDFIP2	PARP9	ITGA1	CEP78	H2AFV	CDCA8	RASGEF1B	ATF3	NUPR1	SPI1	VSIR	LRRC28	COL1A1
43	PASK	BTLA	MAF	CERK	CD3E	AKAP5	IFIH1	IFITM2	ARL4C	CDK1	PLK1	ANXA1	DEDD2	TCEAL9	FCGRT	PASK	BCL11A	PBX1
44	TNFRSF25	NAP1L4	PBXIP1	TLE1	TUBA4A	GOLIM4	SP110	CCL5	BIN2	HIST1H1D	AURKA	MCL1	CD69	MARCKSL1	EGR1	EIF3H	SCAI	GSN
45	FAU	FYB1	F5	EEF1B2	PECAM1	CD27	OAS2	ITGAX	LITAF	HNRNPAB	TROAP	CXCR4	CLK1	SLC40A1	FCGR2A	TXNIP	ATP6AP1L	PDGFRB
46	EEF1D	MIR155HG	SLAMF1	CFH	ARAP2	SNAP47	C19orf66	CD38	TRDC	CENPE	H2AFV	IER5	IER5L	CFH	PLAUR	LDLRAP1	RUFY3	SERPINF1
47	LDLRAP1	PTPN13	BTG3	PERP	ITGA1	RGS1	STAT2	SAMD3	TXK	TK1	KIF2C	MYADM	H3F3B	APOE	CPVL	CMTM8	CD79A	SFRP1
48	CTSL	SESN3	TRAC	PLAT	JAML	HLA-DPA1	LAG3	SLC16A3	MYOM2	ZWINT	MAD2L1	CD8A	NR4A1	GNG11	MS4A7	SCML1	GNA15	LHFPL6
49	ITGA6	BIRC3	IL1R1	KIF5C	CD84	LINC02446	LAP3	CAPN12	GZMM	BIRC5	NUCKS1	PTGER4	CXCR4	CDKN1C	ALDH2	LINC00402	HHIP-AS1	PLAC9
50	PFDN5	TP53INP1	DNPH1	PLCB1	AOAH	SAMSN1	APOL6	CD247	CD300A	PTTG1	RAD21	ZFP36L2	EIF4A2	NBL1	SERPINA1	RIPK2	GSTM3	SERPING1

2 Clusters

2.1 sizes

enframe(sort(table(seu_obj$cluster_label))) %>% 
  mutate(name = ordered(name, levels = rev(name))) %>% 
  ggplot() +
  geom_bar(aes(name, value), stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(y = c("#cells"), x = "cluster")

2.2 UMAP

alpha_lvl <- ifelse(nrow(plot_data) < 20000, 0.2, 0.1)
pt_size <- ifelse(nrow(plot_data) < 20000, 0.2, 0.05)

common_layers_disc <- list(  
  geom_point(size = pt_size, alpha = alpha_lvl),
  NoAxes(),
  guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))),
  labs(color = "")
)

common_layers_cont <- list(  
  geom_point(size = pt_size, alpha = alpha_lvl),
  NoAxes(),
  scale_color_gradientn(colors = viridis(9)),
  guides(color = guide_colorbar())
)

ggplot(plot_data, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  #facet_wrap(~therapy) +
  ggtitle("Sub cluster")

3 Filtering out doublet clusters

3.1 UMAP

my_subtypes <- names(clrs$cluster_label[[coi]])
my_subtypes <- c(my_subtypes, unlist(lapply(paste0("_", 1:3), function(x) paste0(my_subtypes, x)))) %>% .[!str_detect(., "doublet|dissociated")]
my_subtypes <- my_subtypes[my_subtypes %in% unique(seu_obj$cluster_label)]
my_subtypes <- my_subtypes[my_subtypes %in% names(clrs$cluster_label[[coi]])]

if (cell_sort == "CD45+") {
cells_to_keep <- colnames(seu_obj)[seu_obj$cluster_label %in% my_subtypes & !(str_detect(seu_obj$cell_id, "CD45N"))]
}

if (cell_sort == "CD45-") {
cells_to_keep <- colnames(seu_obj)[seu_obj$cluster_label %in% my_subtypes & !(str_detect(seu_obj$cell_id, "CD45P"))]
}

# seu_obj_sub <- subset(seu_obj, cells = cells_to_keep)
# seu_obj_sub <- RunUMAP(seu_obj_sub, dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")
# seu_obj_sub$cluster_label <- seu_obj$cluster_label[colnames(seu_obj) %in% colnames(seu_obj_sub)]
# write_rds(seu_obj_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))
seu_obj_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered.rds"))

plot_data_sub <- as_tibble(FetchData(seu_obj_sub, c(myfeatures, "cluster_label"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cell_id = colnames(seu_obj_sub),
         cluster_label = ordered(cluster_label, levels = my_subtypes),
         )
  
ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label) +
  scale_color_manual(values = clrs$cluster_label[[coi]])

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
  # geom_point(aes(umapharmony_1, umapharmony_2), 
  #            color = "grey90", size = 0.01, 
  #            data = select(plot_data_sub, -tumor_supersite)) +
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

write_tsv(select(plot_data_sub, cell_id, everything(), -umapharmony_1, -umapharmony_2, -contains("RNA_")), paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_embedding.tsv"))

4 sub clustering

big_clusters <- list(
  CD4.T.CXCL13 = "CD4.T.CXCL13", 
  CD4.T.naive = "CD4.T.naive", 
  CD4.T.reg = "CD4.T.reg", 
  CD4.Th17 = "CD4.Th17", 
  CD8.T.CXCL13 = "CD8.T.CXCL13", 
  CD8.T.cytotoxic = "CD8.T.cytotoxic", 
  CD8.T.ISG = "CD8.T.ISG", 
  Cycling.T.NK = "Cycling.T.NK", 
  NK.CD56 = "NK.CD56", 
  NK.cytotoxic = "NK.cytotoxic"
)

cells_to_keep_list <- lapply(big_clusters, function(x) 
colnames(seu_obj_sub)[seu_obj_sub$cluster_label %in% x])

# seu_list <- lapply(cells_to_keep_list, function(x) subset(seu_obj_sub, cells = x))
# 
# preprocess_wrapper <- . %>%
#   FindNeighbors(reduction = "harmony", dims = 1:50) %>%
#   # FindClusters(res = 0.1) %>%
#   # FindClusters(res = 0.2) %>%
#   FindClusters(res = 0.3)
# 
# seu_list <- lapply(seu_list, preprocess_wrapper)
# write_rds(seu_list, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub_clusters.rds"))

seu_list <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub_clusters.rds"))

get_cluster_vector <- function(cluster_label, cluster_res = RNA_snn_res.0.3) {
  cluster_res <- enquo(cluster_res)
  cbind(cell_id = colnames(seu_list[[cluster_label]]), FetchData(seu_list[[cluster_label]], c(as_label(cluster_res)))) %>%
    as_tibble %>%
    mutate(cluster_extended = paste0(cluster_label, "_", !!cluster_res)) %>%
    select(cell_id, cluster_extended) %>% 
    deframe
}

cluster_vector_list <- lapply(big_clusters, get_cluster_vector)

cluster_vector <- unlist(cluster_vector_list, use.names = F) %>% setNames(lapply(cluster_vector_list, names) %>% unlist(use.names = F))

seu_obj$cluster_extended <- cluster_vector[seu_obj$cell_id] %>%
  setNames(seu_obj$cell_id)

seu_obj_sub$cluster_extended <- cluster_vector[seu_obj_sub$cell_id] %>%
  setNames(seu_obj_sub$cell_id)

cluster_extended_uniq <- as.character(na.omit(unique(seu_obj_sub$cluster_extended)))

Idents(seu_obj_sub) <- seu_obj_sub$cluster_extended
Idents(seu_obj) <- seu_obj$cluster_extended

patient_specific_subcluster <- FetchData(seu_obj_sub, c("patient_id", "cluster_extended")) %>% 
  as_tibble() %>% 
  group_by(patient_id, cluster_extended) %>% 
  tally %>% 
  group_by(cluster_extended) %>% 
  mutate(nrel = n/sum(n),
         ntotal = sum(n)) %>% 
  arrange(desc(nrel)) %>% 
  distinct(cluster_extended, .keep_all = T) %>% 
  ungroup

# marker_tbl_extended <- lapply(
#   cluster_extended_uniq,
#   function(x) {
#     marker_tbl <- FindMarkers(seu_obj_sub, ident.1 = x) %>%
#       as_tibble(rownames = "gene") %>%
#       separate(gene, into = c("gene", "drop"), sep = "\\.\\.\\.") %>%
#       select(-drop) %>% 
#       mutate(cluster_extended = x)
#     write_tsv(marker_tbl, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_subcluster_", x, ".tsv"))
#     return(marker_tbl)
# }) %>% bind_rows
#  
# write_tsv(marker_tbl_extended, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated_extended.tsv"))

marker_tbl_extended <- read_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated_extended.tsv"))
 
cluster_label_extended <- c(
  CD4.T.naive_0 = "CD4.T.naive.centr.mem", 
  CD4.T.naive_1 = "CD4.T.effector.memory", 
  CD4.T.CXCL13_0 = "CD4.T.dysfunc.early", 
  CD4.T.CXCL13_1 = "CD4.T.dysfunc.late.1", 
  CD4.T.CXCL13_2 = "CD4.T.dysfunc.late.2", 
  CD4.T.reg_0 = "CD4.T.reg.2", 
  CD4.T.reg_1 = "CD4.T.reg.1", 
  CD4.T.reg_2 = "CD4.T.reg.ISG", 
  CD4.T.reg_3 = "CD4.T.reg.3", 
  CD4.Th17_0 = "CD4.Th17.1", 
  CD4.Th17_1 = "CD4.Th17.2", 
  CD8.T.cytotoxic_0 = "CD8.T.effector.memory", 
  CD8.T.cytotoxic_1 = "CD8.T.cytotoxic", 
  CD8.T.cytotoxic_2 = "dissociated_3", 
  CD8.T.cytotoxic_3 = "CD8.T.naive.centr.mem", 
  CD8.T.cytotoxic_4 = "gd.T.cell", 
  CD8.T.cytotoxic_5 = "doublet.Plasma.cell_1", 
  CD8.T.CXCL13_0 = "CD8.T.dysfunc.early", 
  CD8.T.CXCL13_1 = "CD8.T.dysfunc.late", 
  CD8.T.CXCL13_2 = "NK.reg.CRTAM", 
  CD8.T.ISG_0 = "CD8.T.ISG", 
  CD8.T.ISG_1 = "CD4.T.ISG", 
  CD8.T.ISG_2 = "dissociated_4", 
  NK.CD56_0 = "NK.reg.IGFBP2",
  NK.CD56_1 = "NK.reg.CD56", 
  NK.CD56_2 = "NK.reg.ISG", 
  NK.CD56_3 = "NK.reg.KRT81.KRT86", 
  NK.CD56_4 = "NK.reg.CCL3",
  NK.cytotoxic_0 = "NK.cytotoxic.SPON2.2", 
  NK.cytotoxic_1 = "NK.cytotoxic.GZMH", 
  NK.cytotoxic_2 = "NK.cytotoxic.SPON2.1",
  Cycling.T.NK_0 = "Cycling.CD8.T.1",
  Cycling.T.NK_1 = "Cycling.CD8.T.2", 
  Cycling.T.NK_2 = "Cycling.CD8.T.3", 
  Cycling.T.NK_3 = "Cycling.CD8.T.4", 
  Cycling.T.NK_4 = "Cycling.NK.3", 
  Cycling.T.NK_5 = "Cycling.CD4.T", 
  Cycling.T.NK_6 = "Cycling.NK.2", 
  Cycling.T.NK_7 = "Cycling.NK.1"
)

names(seu_obj$cluster_label) <- colnames(seu_obj)
seu_obj$cluster_label[is.na(seu_obj$cluster_label)] <- "NA"
Idents(seu_obj) <- seu_obj$cluster_label
seu_obj$cluster_label_sub <- cluster_label_extended[seu_obj$cluster_extended]
names(seu_obj$cluster_label_sub) <- colnames(seu_obj)
seu_obj$cluster_label_sub[is.na(seu_obj$cluster_label_sub)] <- "NA"
Idents(seu_obj) <- seu_obj$cluster_label_sub

names(seu_obj_sub$cluster_label) <- colnames(seu_obj_sub)
seu_obj_sub$cluster_label[is.na(seu_obj_sub$cluster_label)] <- "NA"
Idents(seu_obj_sub) <- seu_obj_sub$cluster_label
seu_obj_sub$cluster_label_sub <- cluster_label_extended[seu_obj_sub$cluster_extended]
names(seu_obj_sub$cluster_label_sub) <- colnames(seu_obj_sub)
seu_obj_sub$cluster_label_sub[is.na(seu_obj_sub$cluster_label_sub)] <- "NA"
Idents(seu_obj_sub) <- seu_obj_sub$cluster_label_sub

marker_sheet_extended <- marker_tbl_extended %>%
  mutate(cluster_label_sub = cluster_label_extended[cluster_extended]) %>% 
  group_by(cluster_label_sub) %>%
  mutate(rank = row_number()) %>% 
  slice(1:50) %>% 
  select(rank, gene, cluster_label_sub) %>% 
  arrange(cluster_label_sub) %>% 
  spread(cluster_label_sub, gene)

cluster_n_tbl_full <- seu_obj$cluster_label_sub %>% 
  table() %>% 
  enframe("cluster_label_sub", "cluster_n") %>% 
  mutate(cluster_nrel = cluster_n/sum(cluster_n)) %>% 
  filter(cluster_label_sub != "NA")

marker_tbl_extended_annotated <- marker_tbl_extended %>% 
  mutate(cluster_label_sub = cluster_label_extended[cluster_extended]) %>%
  select(-cluster_extended) %>% 
  left_join(cluster_n_tbl_full, by = "cluster_label_sub") %>% 
  arrange(cluster_label_sub, desc(avg_logFC))

if (cell_sort == "CD45+") {
cells_to_keep <- colnames(seu_obj_sub)[!(str_detect(seu_obj_sub$cluster_label_sub, "dissociated|doublet")) & !(str_detect(seu_obj_sub$cell_id, "CD45N"))]
}

if (cell_sort == "CD45-") {
cells_to_keep <- colnames(seu_obj_sub)[!(str_detect(seu_obj_sub$cluster_label_sub, "dissociated|doublet")) & !(str_detect(seu_obj_sub$cell_id, "CD45P"))]
}


# seu_obj_sub_sub <- subset(seu_obj_sub, cells = cells_to_keep)
# seu_obj_sub_sub <- RunUMAP(seu_obj_sub_sub, dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")
# write_rds(seu_obj_sub_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub.rds"))

seu_obj_sub_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub.rds"))

seu_obj_sub_sub$cluster_label_sub <- setNames(cluster_label_extended[seu_obj_sub_sub$cluster_extended], seu_obj_sub_sub$cell_id)

write_rds(seu_obj_sub_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_sub.rds"))

formattable::formattable(marker_sheet_extended)

rank	CD4.T.dysfunc.early	CD4.T.dysfunc.late.1	CD4.T.dysfunc.late.2	CD4.T.effector.memory	CD4.T.ISG	CD4.T.naive.centr.mem	CD4.T.reg.1	CD4.T.reg.2	CD4.T.reg.3	CD4.T.reg.ISG	CD4.Th17.1	CD4.Th17.2	CD8.T.cytotoxic	CD8.T.dysfunc.early	CD8.T.dysfunc.late	CD8.T.effector.memory	CD8.T.ISG	CD8.T.naive.centr.mem	Cycling.CD4.T	Cycling.CD8.T.1	Cycling.CD8.T.2	Cycling.CD8.T.3	Cycling.CD8.T.4	Cycling.NK.1	Cycling.NK.2	Cycling.NK.3	dissociated_3	dissociated_4	doublet.Plasma.cell_1	gd.T.cell	NK.cytotoxic.GZMH	NK.cytotoxic.SPON2.1	NK.cytotoxic.SPON2.2	NK.reg.CCL3	NK.reg.CD56	NK.reg.CRTAM	NK.reg.IGFBP2	NK.reg.ISG	NK.reg.KRT81.KRT86
1	CXCL13	CXCL13	CXCL13	IL7R	IFIT3	CCR7	RTKN2	TNFRSF4	TNFRSF4	ISG15	KLRB1	LST1	GZMH	CXCL13	CXCL13	GZMK	IFIT3	CD8B	STMN1	STMN1	STMN1	CXCL13	CXCL13	SPON2	CENPF	CENPF	GZMK	XIST	IGHG3	TRGC2	FGFBP2	SPON2	FGFBP2	CCL3	GNLY	XCL2	TYROBP	ISG15	GNLY
2	MAF	NMB	PDCD1	CD40LG	ISG15	RPS8	FOXP3	IL2RA	TNFRSF18	SAT1	IL7R	IL4I1	CD8A	CCL4L2	KRT86	CD8A	ISG15	IL7R	MKI67	TUBA1B	MKI67	STMN1	HLA-DRA	FCGR3A	ASPM	MKI67	MALAT1	MALAT1	IGHG1	XCL1	GZMH	FGFBP2	FCGR3A	AREG	AREG	GZMB	KLRC1	GNLY	KLRC1
3	IL6ST	CPM	CTLA4	ANXA1	MX2	RPS6	LTB	FOXP3	TNFRSF9	HERC5	IL4I1	TNFSF13B	CD8B	HAVCR2	TNFRSF9	CD8B	IFIT1	CD8A	TUBA1B	HIST1H4C	TOP2A	TYMS	HLA-DRB1	IGFBP7	MKI67	STMN1	PTPRC	NEAT1	IGLC3	ZNF683	CX3CR1	FCGR3A	SPON2	CCL4L2	KLRC1	CRTAM	FCER1G	IFIT3	KRT81
4	NR3C1	NR3C1	MAF	GPR183	IFIT2	RPS3A	UGP2	BATF	IL2RA	MX1	LTB	PCDH9	CCL5	MIR155HG	GZMB	KLRG1	MX1	YBX3	HMGB2	MKI67	ASPM	TUBA1B	GZMB	KLRF1	TOP2A	TOP2A	EDF1	MT-ND3	IGLC2	KLRC2	NKG7	CCL3	KLRF1	XCL2	ITGAX	XCL1	TRDC	FCER1G	FCER1G
5	FKBP5	IGFL2	ZBED2	KLF2	IFI44L	SELL	TBC1D4	TNFRSF18	BATF	TNFRSF4	CCR6	SPINK2	HLA-DRB1	PTMS	RBPJ	ITM2C	IFIT2	ANXA1	TOP2A	TUBB	TYMS	PCNA	IFNG	CX3CR1	KNL1	ASPM	RPS29	MT-ND1	IGHG4	LINC02446	TRGC2	PRF1	PRF1	XCL1	FCER1G	FABP5	GNLY	MX1	AREG
6	CD40LG	FKBP5	RBPJ	S100A4	IFIT1	EEF1B2	IL2RA	LTB	CTLA4	SPATS2L	TNFSF13B	KIT	GZMA	CD8A	SRGAP3	CCL4	RSAD2	RPS12	TUBB	TYMS	HELLS	PCLAF	FABP5	AKR1C3	NUSAP1	TUBA1B	RPL36AL	MT-CO3	IGKC	HOPX	KLF2	KLRF1	GNLY	FOS	XCL1	TNFRSF9	AREG	IFIT1	TNFRSF18
7	CD4	ITM2A	DUSP4	FOS	RSAD2	RPS12	PMAIP1	CTLA4	GADD45A	IFI6	LST1	SCN1B	GZMK	LAG3	FAM3C	DUSP2	IFI6	TOB1	CENPF	TOP2A	CLSPN	DUT	COTL1	GINS2	CENPE	HMGB2	RPL13A	MT-ND2	JCHAIN	XCL2	KLRG1	PTGDS	KLRD1	IER2	TYROBP	HSP90AB1	IGFBP2	TYROBP	KRT86
8	AC004585.1	ICA1	CD4	LTB	MX1	RPL34	IKZF2	GADD45A	TNFRSF1B	IL2RA	AQP3	LTC4S	HLA-DPB1	GZMB	TIGIT	PPP1R14B	ISG20	FOS	HMGN2	HMGB2	PCNA	FABP5	GAPDH	FGFBP2	TPX2	GNLY	ATP5MC2	MT-CYB	IGHGP	CTSW	PLEK	CX3CR1	NKG7	CCL4	TRDC	PKM	KRT81	IFI27	TRDC
9	CORO1B	IL6ST	CD40LG	TPT1	IL7R	RPL32	IL32	RTKN2	ICOS	CTLA4	CCL20	IL1R1	PTMS	IFNG	SPRY1	TRGC2	MX2	VIM	PTTG1	NUSAP1	UBE2C	CLSPN	NASP	MLC1	DLGAP5	TUBB	CALM1	MT-ATP6	IGHG2	CD63	PRF1	GZMB	CLIC3	EGR1	XCL2	GAPDH	KLRD1	IFITM3	TYROBP
10	ITM2A	TSHZ2	MIR155HG	NFKBIA	IFI6	KLF2	CTLA4	SOX4	PIM3	OAS1	CTSH	LINC00299	ZNF683	CCL3	IFNG	CST7	SAMD9L	JUNB	IL2RA	CENPF	PCLAF	HELLS	LAG3	S1PR5	GTSE1	NUSAP1	SRP14	MT-CO1	CD8A	CD7	ZEB2	IGFBP7	TYROBP	NFKBIA	KLRD1	NME1	CLIC3	KLRC1	XCL1
11	TSHZ2	TNFRSF4	PRDM1	TIMP1	TNFSF10	RPL9	TIGIT	TNFRSF1B	FOXP3	IFIT3	RORA	GSN	HLA-DPA1	CXCR6	MIR155HG	DTHD1	IFI44L	GZMK	PCLAF	PCLAF	GTSE1	NASP	MCM5	TYMS	KIF14	UBE2C	COMMD6	PTPRC	MZB1	TRGC1	PRSS23	PLAC8	GZMB	CCL3L1	CMC1	MIR155HG	KLRB1	AREG	CTSW
12	PDCD1	AHI1	TNFRSF18	S100A10	STAT1	IL7R	ARID5B	TIGIT	IL1R2	IFIT1	TNFRSF25	PLAT	CXCR6	HLA-DRB1	CD8A	CCL4L2	HERC5	CCR7	NUSAP1	HIST1H1B	ORC6	MCM3	PTMS	E2F1	HJURP	IGFBP2	RPS25	MT-CO2	CCL5	NCR3	KLRD1	GNLY	PLAC8	DUSP2	CTSW	SRM	CEBPD	IFI44L	ZNF683
13	TOX2	AC004585.1	ITM2A	JUN	HERC5	RPS13	SAT1	TBC1D4	LTB	FOXP3	SLC4A10	CD300LF	HLA-DRA	CCL5	LAG3	CRTAM	OAS1	ZFP36L2	ASPM	HMGN2	CDT1	MKI67	MIR155HG	MCM7	POLQ	TPX2	PFDN5	H3F3B	IGHA1	CD52	FCGR3A	FCER1G	CX3CR1	BTG2	NCAM1	NPW	KRT86	RSAD2	XCL2
14	CH25H	LIMS1	RGS1	S100A11	ISG20	RPL21	BATF	LINC01943	TBC1D4	TYMP	DPP4	AFF3	ITM2C	RBPJ	PHLDA1	GZMM	MT2A	S100A6	H2AFZ	ASPM	BRCA1	MCM5	HAVCR2	MYOM2	BUB1B	KLRC1	RPS15A	CCNI	CCL4	TRDC	MYBL1	NKG7	FCER1G	NR4A1	HOPX	KDM6B	TXK	CD38	HOPX
15	PTPN13	BIRC3	LINC01943	VIM	OAS1	TCF7	CORO1B	SAT1	PKM	LTB	MYBL1	DLL1	COTL1	GZMH	ENTPD1	CMC1	STAT1	MT-ND1	TYMS	H2AFZ	ZWINT	MCM7	CD8A	MCM4	HMMR	HMGN2	HINT1	GMFG	GZMK	KLRD1	ADGRG1	AKR1C3	PTGDS	FOSB	LINC00996	TSPAN17	XCL1	IFIT2	RGS16
16	CTSW	PGM2L1	TNFRSF4	CTSL	OAS3	MAL	MIR4435-2HG	TNFRSF9	BIRC3	RTKN2	S100A4	LIF	HLA-DQA1	HLA-DRA	HAVCR2	LYAR	TNFSF10	KLF2	CEP55	RRM2	ATAD2	MCM4	HELLS	CHEK1	KIF11	SMC4	RPS27	COX5B	CD8B	KIR3DL2	EFHD2	AREG	PLEK	CD69	KLRF1	ZBED2	IL2RB	IFI6	CD7
17	CCL5	CD200	ADAM19	RPL34	EIF2AK2	RPL39	TNFRSF4	UGP2	REL	EPSTI1	CD40LG	IL23R	IFNG	PHLDA1	CCL3	CCL5	OASL	LYAR	DUT	SMC4	RAD51AP1	TUBB	MCM6	CLSPN	CIT	FCER1G	FTH1	CALM2	RGS1	TRG-AS1	FLNA	CCL4	EFHD2	TYROBP	IFITM3	GOLIM4	MATK	TRDC	LGALS1
18	NKG7	SESN3	CCDC50	EEF1A1	IFI44	RPL22	AC005224.3	IKZF2	ENTPD1	STAT1	TPT1	RORC	CCL4	PDCD1	CTLA4	EOMES	PLSCR1	CXCR4	PCNA	PCNA	NUSAP1	GAPDH	ACTB	STMN1	SGO2	TYMS	HSPA8	PSMB1	GZMA	CCL5	GNLY	KLRD1	CST7	NR4A2	TXK	PAICS	NKG7	TNFSF10	CD63
19	KLRD1	RNF19A	TSHZ2	RPLP0	GBP1	RPL7	PBXIP1	PMAIP1	TYMP	ICOS	RPLP1	TGM2	APOBEC3G	HLA-DPA1	JAML	SH2D1A	GBP1	NELL2	TNFRSF4	DUT	DLGAP5	GINS2	MCM4	CHST2	C21orf58	GSTP1	DBI	TPM3	HLA-DPA1	TMSB4X	CST7	PLEK	CCL3	GADD45B	IL2RB	MRTO4	XCL2	ISG20	KLRC2
20	ST8SIA1	CHN1	CD2	RPS18	CCR7	RPL5	STAM	IL32	SAT1	ISG20	EEF1A1	TNFSF11	SH3BGRL3	CD8B	CCND2	YBX3	SAMD9	S100A10	TPX2	ATAD2	CDK1	HMGB2	STMN1	CTBP2	KIF23	KRT81	FAU	SOD1	CXCR6	LINC01871	HOPX	CLIC3	IGFBP7	KLRD1	KLRC2	GEM	CD63	GZMB	TMIGD2
21	LIMS1	GNG4	ETV7	RPS12	EPSTI1	RPS18	GADD45A	ICOS	IL1R1	MX2	B3GALT2	ZG16B	CD52	HLA-DPB1	PTMS	GZMA	EIF2AK2	MT-ND2	FOXP3	CLSPN	BIRC5	RANBP1	PCNA	FCER1G	AC007240.1	TRDC	PARK7	ATP5MPL	IKZF3	GNLY	C12orf75	EFHD2	ADGRG1	TRDC	MATK	EGR2	TMIGD2	SAMD9L	SRGAP3
22	SMCO4	MAF	MYO7A	RPL39	MT2A	RPL37	BTG3	LAYN	SDC4	TIGIT	RORC	ENPP1	CD3D	VCAM1	LINC01871	TC2N	LAG3	TNFRSF18	ENO1	UBE2C	TUBA1B	HLA-DRA	CLSPN	PLAC8	KIF15	H2AFZ	RPLP2	AES	LTB	ZFP36L2	GZMB	MYOM2	HOPX	KLRC1	SELL	CD72	CMC1	BST2	MATK
23	CD8A	BATF	TOX2	JUNB	SAMD9L	RPL13	BIRC3	LINC02099	LINC01943	RSAD2	LTK	PTGER3	LAG3	LINC01871	GOLIM4	CXCR4	IFI35	LMNA	BATF	HMGB1	NASP	MCM6	GINS2	RAMP1	NCAPG	TUBB4B	C9orf16	GNG5	APOC1	HCST	C1orf21	TYROBP	CD247	FCER1G	CD7	METTL1	CTSW	CLIC3	CXXC5
24	CTLA4	TOX2	CD200	TOB1	XAF1	RPS5	HPGD	CTSC	NAMPT	TNFRSF1B	CFH	IRAK3	C12orf75	JAML	GAPDH	PPP2R5C	EPSTI1	TSC22D3	CTLA4	HIST1H1D	MCM4	DEK	CDCA7	ORC6	BRCA2	AREG	RPL23A	SAP18	RPL22	IKZF2	LITAF	ADGRG1	AREG	NFKBID	GZMB	GFOD1	CD7	PLSCR1	GSTP1
25	CCL4	CH25H	BHLHE40-AS1	AQP3	TMEM123	RPL3	SELL	DUSP4	LTA	UGP2	NCR3	CA2	CD74	TNIP3	TNS3	A2M-AS1	OAS3	FOSB	GAPDH	TPX2	GINS2	H2AFZ	CHEK1	MCM2	ECT2	CENPE	RPL28	CIB1	PABPC1	FCRL6	PLAC8	CST7	ZEB2	ZFP36	SRGAP3	RGS16	CCL3	OAS1	IL2RB
26	FYB1	THADA	PTPN13	RPL11	DDX58	EEF1A1	ICOS	CD27	DNPH1	IFI44L	FAM241A	KIAA1211L	CD3G	TNFSF4	LAYN	GZMH	XAF1	ZFP36	HIST1H1B	H2AFX	CENPU	HMGN2	MCM2	TTC38	CKAP2L	UBE2S	COX7C	COPE	HLA-DPB1	IL32	S100A4	HOPX	AKR1C3	CEBPD	CEBPD	PARK7	LAT2	IFI35	LAT2
27	RUNX3	CD84	PDE7B	TNFAIP3	CMPK2	RPL11	AC133644.2	ARID5B	TIGIT	SAMD9	CEBPD	SLC4A10	THEMIS	CD63	NBL1	PDCD4	USP18	PPP1R14B	SMC4	SMC2	FEN1	SMC2	TYMS	GZMB	ARHGAP11B	CKS1B	CHCHD2	TRMT112	TPT1	NUCB2	SPON2	CD247	KLRB1	MAFF	CD63	NDFIP2	GSTP1	KLRD1	CSF1
28	RNF19A	BTLA	NKG7	TNFRSF25	CD40LG	RPS23	NABP1	CORO1B	TFRC	TNFRSF18	RPLP0	COL4A4	LINC02446	HLA-DQA1	MYO1E	HLA-DPB1	LY6E	DUSP1	SMC2	TMPO	TUBB	DNMT1	CDT1	PCNA	ANLN	TNFRSF18	RPL4	PSMB9	CD74	CD9	FCRL6	PRSS23	PTPN12	IGFBP2	NKG7	IFNG	HOPX	EPSTI1	SPRY1
29	PRF1	CORO1B	IFNG	FXYD5	SAMD9	RPL23	CD27	ENTPD1	F5	SAMD9L	RPS12	NRP1	HLA-DQB1	ID2	RGS2	F2R	CMPK2	CD55	BIRC5	CKS1B	CHEK1	TPI1	MALAT1	PLEK	ESPL1	HIST1H1B	RPS24	BRK1	GZMH	CD8A	ADRB2	ABHD17A	PRSS23	KDM6B	TNFRSF18	COL6A3	CCL5	IRF7	CAPG
30	RGS1	METTL8	NR3C1	RPL3	LY6E	TPT1	IKZF4	MAGEH1	PMAIP1	TBC1D4	ERN1	B3GALT5	LGALS1	FABP5	MYO7A	SLAMF7	NT5C3A	CD69	PKM	UBE2S	KNL1	ATAD2	HLA-DQA1	MCM5	NUF2	KPNA2	COX8A	RBX1	IGHM	KLRC4	ARL4C	CTSW	C1orf21	IER3	IFITM2	CXCL13	KLRF1	CEBPD	ITGA1
31	GZMB	CCDC50	LIMS1	RPL9	GPR183	RPLP0	PHTF2	S100A4	CCND2	BATF	KIF5C	PDZK1	HLA-DRB5	LINC02446	AKAP5	TNFSF9	IFI44	GPR183	RRM2	KNL1	SMC2	BRCA1	MYO7A	HES6	KIF4A	DLGAP5	RPS16	HSPA8	RPL8	MATK	S1PR5	ZEB2	TXK	CD160	CD300A	TNFRSF18	SH2D1B	XCL2	CLIC3
32	IL6R	SMCO4	RNF19A	RPL32	USP18	RPL30	PHACTR2	BIRC3	DUSP4	USP18	RPL13	EMID1	S100A6	TNFRSF9	DGKH	STK17A	HELZ2	TYROBP	PMAIP1	TUBB4B	MCM5	HNRNPAB	MKI67	PRF1	KIFC1	CKS2	RPL35A	S100A11	CD3E	CRTAM	KLF3	CHST2	S1PR5	REL	GSTP1	RGCC	PRF1	HERC5	LDLRAD4
33	CST7	CD4	RDH10	RPS3A	DDX60L	RPS28	RGS1	CD4	LINC02099	COX5A	SPOCK2	HOXA5	FABP5	TIGIT	SNX9	PKM	IRF7	FCER1G	TK1	HIST1H1C	KIF11	SLBP	HLA-DPA1	FEN1	BUB1	XCL2	RPL7	VAMP8	HLA-DQB1	S100A6	PATL2	PTPN12	CEBPD	GNLY	KRT81	CD82	ITGA1	NT5C3A	KIR2DL4
34	SLC9A9	ARMH1	GAPDH	RPS6	TRIM22	LEF1	DUSP16	MAF	CCL20	SOX4	JAML	LINGO4	ANXA5	ITGAE	DUSP4	COX5A	CD38	CCL3	JPT1	EZH2	TPX2	SNRPB	SH3BGRL3	CDT1	KIF20B	SMC2	RPS7	POMP	WNK1	CAPG	ANXA1	MYBL1	TRDC	KLRF1	SPTSSB	MALAT1	SRGAP3	IL2RB	IGFBP2
35	BTLA	PHACTR2	SNX9	CD4	HELZ2	RPS20	TRAC	SPOCK2	SH2D2A	IKZF2	TIGIT	NECTIN2	CCL4L2	LYST	LINC01480	EIF2AK2	TRIM22	CTSC	ANP32B	DEK	HIST1H1B	GZMB	MCM3	ADGRG1	FBXO43	TMPO	RPS11	EDF1	CD4	IFITM2	TGFBR3	TTC38	CTSW	CLIC3	IGFBP2	RBPJ	MCTP2	LY6E	TXK
36	TP53INP1	TNFRSF25	IL6ST	RPL36	ANXA1	RPL38	CARD16	IL1R1	ZBTB32	PMAIP1	RPL34	HOXA7	TMSB4X	APOBEC3C	ITGAE	ACP5	DDX58	RGCC	TMPO	CDK1	GMNN	COTL1	PTTG1	CCDC28B	NCAPG2	CDK1	ATP5MG	NDUFS5	HLA-DRA	MAF	CTSW	S1PR5	GZMH	TXK	CAPN12	RANBP1	MAPK1	XAF1	LINC01871
37	TNFSF8	ARID5B	CORO1B	RPL30	NT5C3A	RPS16	LINC01943	GLRX	ID3	CTSC	S100A6	PRAM1	LGALS3	CD27	RGS1	ADGRG1	PARP14	AC020916.1	CENPE	BIRC5	MCM7	DNAJC9	MCM7	GNLY	PRC1	CTSW	MYL12B	UQCR11	RPS13	CD4	S100A6	ITGB2	TTC38	MATK	CXXC5	EIF5A	CXXC5	XCL1	CD9
38	GZMH	CTLA4	ITGA2	RPS8	PLSCR1	RPS4X	S100A4	CARD16	SOX4	NT5C3A	RPL32	B4GALNT1	ITGA1	HLA-DMA	CD63	IFI35	CD8A	CTLA4	UBE2C	H2AFV	RRM1	HLA-DRB1	RBPJ	ASCL2	DIAPH3	KNL1	RPL39	LDHB	CD3G	GLUL	TTC16	XBP1	ABHD17A	KLRB1	IL18	DCTPP1	KLRC2	PRF1	SLC16A3
39	GNLY	ZBED2	SMCO4	RPS27A	IFIH1	RPL36	TNFRSF1B	TYMP	PTP4A3	TNFRSF9	HLA-DPB1	LPAR1	TRGC2	CRTAM	SNAP47	CSF1	DDX60	CD4	CLSPN	TK1	WDR76	CENPX	CD8B	PTPN12	SMC4	KRT86	UBB	HNRNPA1	RPL29	PRMT9	A2M-AS1	CD160	PTGDR	NKG7	ID2	ENO1	SAMD3	SAMD9	CCL5
40	METTL8	TNFRSF18	SAMSN1	RPSA	IFI35	RPLP2	MAGEH1	MIR4435-2HG	MIR155HG	IFI44	TLE1	SLCO2A1	APOBEC3C	NKG7	CLNK	NDFIP2	PPM1K	DNAJB1	CKS2	CENPM	RRM2	IFNG	PDCD1	PCLAF	RPL10	PCLAF	RPS21	ANP32B	TNFRSF18	SCML4	TTC38	LAIR2	MYBL1	EGR2	PTGDR	TPI1	LINC00996	MX2	NCR3
41	ADAM19	RILPL2	SRGN	RPL13	IRF7	RPS27	CD4	DNPH1	ARID5B	ARID5B	TMIGD2	CD33	SIT1	CTLA4	TNIP3	LGALS3	IFIH1	MYBL1	ACTB	ZWINT	CENPF	CHEK1	BRCA1	CDC45	RPS27	TYROBP	RPL13	ABRACL	LARP7	LYAR	ASCL2	KLRB1	CMC1	IL2RB	KLRB1	CAMK1	CD247	OASL	PRMT9
42	CD2	AGFG1	ICOS	EEF1B2	DDX60	RPS2	SPOCK2	PHACTR2	CD27	LINC01943	FKBP11	APOL4	JAML	RGS1	ADAMTS6	TMEM123	RNF213	CD38	CALM3	HIST1H1E	UBE2T	CCL3	ITGA2	GMNN	BRCA1	XCL1	ANP32B	TRAPPC1	DUSP23	CTSA	VCL	TXK	CHST2	MAP3K8	SH2D1B	TNIP3	PTGDR	USP18	CLNK
43	SPOCK2	CYSLTR1	KLRD1	RPL38	OAS2	RPL12	CCNG2	IL1R2	ENO1	GBP5	HLA-DRB1	RASSF8	LINC01871	ACP5	BCL2L11	IL6ST	DDX60L	TNFRSF4	ATAD2	FAM111B	CDCA7	CDT1	CD63	CENPU	RPL41	HMMR	COX6C	ATP5F1D	DUSP2	IFITM3	LILRB1	PTGDR	CCL4	IRF8	C1orf162	PIM3	NCAM1	EIF2AK2	CTSA
44	GZMA	CD40LG	CD109	RPLP1	PNPT1	SNHG8	FCMR	PHLDA1	IKZF2	OAS3	RPS27A	TLE1	CD2	SNAP47	NDFIP2	CD247	OAS2	TNFAIP3	ANP32E	KIF11	CDKN3	RAN	CCL5	C1orf21	RPL30	TMIGD2	RPS28	SNRPB	SELL	PTGDR	FGR	CD300A	SH2D1B	KRT81	PLAC8	NAMPT	STARD3NL	ERICH3	KLRB1
45	CD8B	BHLHE40-AS1	CD3D	RPL5	PPM1K	RPS27A	TLK1	STAM	CD4	CD4	KLRD1	AL034397.3	CST7	ITM2A	RAB27A	TPI1	PARP9	HLA-DRA	ICOS	FABP5	SGO2	PKM	CCL3	FGR	RPL13	CDKN3	RPL36	SEPT7	KLF2	SPRY1	PLEKHG3	FGR	CEP78	SRGN	KIR2DL4	TNFRSF1B	CAPN12	CCL3	CTSD
46	MATK	IGFBP4	SLAMF1	RPS16	PARP9	RPL19	SMAP2	SLAMF1	LAIR2	RNF213	RUNX2	STAC	CCR5	SAMSN1	TNFRSF1B	ANKRD28	PNPT1	JUN	CKS1B	DLGAP5	USP1	ENO1	HLA-DMA	LYN	RPS27A	H2AFX	LDHB	ZFAS1	CXCR4	CTLA4	COL6A2	GSAP	XBP1	LAT2	MCTP2	ZNF593	TNFRSF18	IFI44	NDFIP2
47	CD84	ELMO1	RPL23A	CD52	LGALS3BP	PABPC1	LAYN	LAIR2	SYNGR2	SPOCK2	RPL10	BCAS1	CRTAM	CCND2	ID2	CD226	GZMK	NCF1	PHF19	GTSE1	FANCI	SNRPD1	GZMH	KIR2DL3	RPS12	CD63	SUMO2	UQCRH	DNAJC7	DDAH2	CES1	CTBP2	SAMD3	ZBTB16	ZNF683	HSPA5	AOAH	CD7	ADGRG3
48	ARID5B	MS4A6A	AHI1	LMNA	LTB	RPS14	PIM2	F5	SERPINB9	ENTPD1	HLA-DPA1	CSF2	GRAP2	APOBEC3G	LYST	TMEM173	BST2	TRDC	TNFRSF18	TMEM106C	MT-ND4	NUDT1	APOBEC3G	ITGAM	PRR11	CCNB1	RPL14	SEC61G	TRBC2	GUK1	LINC02384	GK5	MATK	GRASP	MMP25-AS1	FAM3C	GZMA	OAS3	CD247
49	CD63	ST8SIA1	ARID5B	RPL37	LAMP3	RPS21	SAMHD1	TRAC	NR4A1	CD27	RPS13	TOX2	IDH2	CD74	AHI1	TPM4	SP110	HAVCR2	UBE2S	ANP32B	ATAD5	PGAM1	ZNF683	MCM10	RPS14	GTSE1	RPS19	PPP1CA	SERPINB6	CD28	COTL1	C1orf21	ITGB2	NR4A3	MAFF	SLC7A5	APOBEC3G	KLRC2	GZMB
50	HOPX	MAGEH1	TYMP	RPS29	GZMA	RPL35A	ICA1	PHTF2	VDR	LY6E	GZMB	ALDOC	GZMM	DUSP4	ZEB2	USP18	C19orf66	MYADM	RANBP1	HNRNPAB	CENPE	ATAD5	ANXA5	CD300A	TPT1	BIRC5	ATP5F1E	RHOA	TYROBP	CCR7	CD27	MTSS1	CD300A	EGR3	XBP1	PGAM1	METRNL	WARS	CKLF

write_tsv(marker_sheet_extended, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_sheet_full.tsv"))

write_tsv(marker_sheet_extended, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_sheet_full.tsv"))

write_tsv(marker_tbl_extended_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_marker_table_annotated_full.tsv"))

write_tsv(marker_tbl_extended_annotated, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/supplementary_tables/", coi, "_marker_table_annotated_full.tsv"))

plot_data_sub_sub <- as_tibble(FetchData(seu_obj_sub_sub, c(myfeatures, "cluster_label_sub"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(cell_id = colnames(seu_obj_sub_sub),
         cluster_label_sub = ordered(cluster_label_sub, levels = names(clrs$cluster_label_sub[[coi]])),
         ) %>%
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite)))
  
if (cell_sort == "CD45+") {
plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
}

if (cell_sort == "CD45-") {
plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
}

ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label_sub)) + 
  common_layers_disc +
  # facet_wrap(~cluster_label_sub) +
  scale_color_manual(values = clrs$cluster_label_sub[[coi]]) +
  ggtitle("Sub cluster") 

4.1 add module scores and pathway scores

signature_modules_cd8 <- read_excel("_data/small/signatures/SPECTRUM v5 sub cluster markers.xlsx", sheet = 2, skip = 1, range = "M2:P100") %>% 
  gather(module, gene) %>% 
  na.omit() %>% 
  group_by(module) %>% 
  do(gene = c(.$gene)) %>% 
  {setNames(.$gene, .$module)} %>% 
  setNames(paste0(names(.), ".module"))

signature_modules_cin <- yaml::read_yaml("_data/small/signatures/cgas_sting.yaml") %>% 
  .[lapply(., length) > 1] %>% 
  setNames(paste0(names(.), ".module"))

## compute expression module scores
for (i in 1:length(signature_modules_cd8)) {
  seu_obj_sub_sub <- AddModuleScore(seu_obj_sub_sub, features = signature_modules_cd8[i], name = names(signature_modules_cd8)[i])
  seu_obj_sub_sub[[names(signature_modules_cd8)[i]]] <- seu_obj_sub_sub[[paste0(names(signature_modules_cd8)[i], "1")]]
  seu_obj_sub_sub[[paste0(names(signature_modules_cd8)[i], "1")]] <- NULL
  print(paste(names(signature_modules_cd8)[i], "DONE"))
}

## [1] "CD8.Cytotoxic.module DONE"
## [1] "CD8.Dysfunctional.module DONE"
## [1] "CD8.Naive.module DONE"
## [1] "CD8.Predysfunctional.module DONE"

for (i in 1:length(signature_modules_cin)) {
  seu_obj_sub_sub <- AddModuleScore(seu_obj_sub_sub, features = signature_modules_cin[i], name = names(signature_modules_cin)[i])
  seu_obj_sub_sub[[names(signature_modules_cin)[i]]] <- seu_obj_sub_sub[[paste0(names(signature_modules_cin)[i], "1")]]
  seu_obj_sub_sub[[paste0(names(signature_modules_cin)[i], "1")]] <- NULL
  print(paste(names(signature_modules_cin)[i], "DONE"))
}

## [1] "CIN.module DONE"
## [1] "CIN.responsive.Noncanonical.NFKB.targets.up.module DONE"
## [1] "CIN.responsive.Noncanonical.NFKB.targets.down.module DONE"
## [1] "Noncanonical.NFKB.regulators.up.module DONE"
## [1] "Noncanonical.NFKB.regulators.down.module DONE"
## [1] "Canonical.NFKB.regulators.module DONE"
## [1] "Interferon.regulators.module DONE"
## [1] "EMT.module DONE"
## [1] "Inflammation.module DONE"
## [1] "Motility.and.migration.module DONE"
## [1] "ISG.module DONE"
## [1] "SASP.module DONE"
## [1] "Myofibroblast.module DONE"
## [1] "Fibroblast.module DONE"
## [1] "Mesenchymal.stem.cell.module DONE"
## [1] "Pericytes.module DONE"
## [1] "Smooth.muscle.cell.module DONE"
## [1] "ER.stress.module DONE"
## [1] "IFNg.signaling.module DONE"

## compute progeny scores
progeny_list <- seu_obj_sub_sub@assays$RNA@data[VariableFeatures(seu_obj_sub_sub),] %>%
  as.matrix %>%
  progeny %>%
  as.data.frame %>%
  as.list

names(progeny_list) <- make.names(paste0(names(progeny_list), ".pathway"))

for (i in 1:length(progeny_list)) {
  seu_obj_sub_sub <- AddMetaData(seu_obj_sub_sub,
                                 metadata = progeny_list[[i]],
                                 col.name = names(progeny_list)[i])
}

write_rds(seu_obj_sub_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_annotated.rds"))

seu_obj_sub_sub <- read_rds(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_annotated.rds"))

seu_obj_sub_sub$cluster_label_sub <- setNames(cluster_label_extended[seu_obj_sub_sub$cluster_extended], seu_obj_sub_sub$cell_id)

write_rds(seu_obj_sub_sub, paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_processed_filtered_annotated.rds"))

4.2 marker heatmap

marker_top_tbl <- marker_sheet_extended[,-1] %>% 
  slice(1:10) %>% 
  as.list %>% 
  .[!str_detect(names(.), "doublet|dissociated")] %>% 
  enframe("cluster_label_x", "gene") %>% 
  unnest(gene)

plot_data_markers <- as_tibble(FetchData(seu_obj_sub_sub, c("cluster_label_sub", myfeatures, unique(marker_top_tbl$gene)))) %>% 
  gather(gene, value, -c(1:(length(myfeatures)+1))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cluster_label_sub = ordered(cluster_label_sub, levels = names(clrs$cluster_label_sub[[coi]]))) %>% 
  group_by(cluster_label_sub, gene) %>% 
  summarise(value = mean(value, na.rm = T)) %>% 
  group_by(gene) %>% 
  mutate(value = scales::rescale(value)) %>% 
  left_join(marker_top_tbl, by = "gene") %>% 
  mutate(cluster_label_x = ordered(cluster_label_x, levels = rev(names(clrs$cluster_label_sub[[coi]]))))

ggplot(plot_data_markers) +
  geom_tile(aes(gene, cluster_label_sub, fill = value)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  facet_wrap(~cluster_label_x, scales = "free") +
  scale_fill_gradientn(colors = viridis(9)) +
  labs(fill = "Scaled\nexpression") +
  theme(aspect.ratio = 1,
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank())

# ggsave(paste0("_fig/002_marker_heatmap_", coi, ".pdf"), width = nrow(marker_top_tbl)/6, height = 5)

4.3 composition

4.3.1 per site

comp_site_tbl <- plot_data_sub_sub %>%
  filter(!is.na(tumor_supersite)) %>% 
  group_by(cluster_label_sub, tumor_supersite) %>%
  tally %>%
  group_by(tumor_supersite) %>%
  mutate(nrel = n/sum(n)*100) %>%
  ungroup

pnrel_site <- ggplot(comp_site_tbl) +
  geom_bar(aes(tumor_supersite, nrel, fill = cluster_label_sub),
           stat = "identity", position = position_stack()) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(fill = "Cluster", y = "Fraction [%]", x = "") +
  scale_fill_manual(values = clrs$cluster_label_sub[[coi]])

pnabs_site <- ggplot(comp_site_tbl) +
  geom_bar(aes(tumor_supersite, n, fill = cluster_label_sub),
           stat = "identity", position = position_stack()) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  labs(fill = "Cluster", y = "# cells", x = "") +
  scale_fill_manual(values = clrs$cluster_label_sub[[coi]])

plot_grid(pnabs_site, pnrel_site, ncol = 2, align = "h")

# ggsave(paste0("_fig/02_deep_dive_", coi, "_comp_site.pdf"), width = 8, height = 4)

4.3.2 per sample

comp_tbl_sample_sort <- plot_data_sub_sub %>% 
  group_by(tumor_subsite, tumor_supersite, patient_id, consensus_signature, therapy, cluster_label_sub) %>% 
  tally %>% 
  group_by(tumor_subsite, tumor_supersite, patient_id, consensus_signature, therapy) %>% 
  mutate(nrel = n/sum(n)*100,
         nsum = sum(n),
         log10n = log10(n),
         label_supersite = "Site",
         label_mutsig = "Signature",
         label_therapy = "Rx") %>% 
  ungroup %>% 
  arrange(desc(therapy), tumor_supersite) %>% 
  mutate(tumor_subsite_rx = paste0(tumor_subsite, "_", therapy)) %>% 
  mutate(tumor_subsite = ordered(tumor_subsite, levels = unique(tumor_subsite)),
         tumor_subsite_rx = ordered(tumor_subsite_rx, levels = unique(tumor_subsite_rx))) %>% 
  arrange(patient_id) %>% 
  mutate(label_patient_id = ifelse(as.logical(as.numeric(fct_inorder(as.character(patient_id)))%%2), "Patient1", "Patient2"))

sample_id_x_tbl <- plot_data_sub %>% 
  mutate(sort_short_x = cell_sort) %>% 
  distinct(patient_id, sort_short_x, tumor_subsite, therapy, sample) %>% 
  unite(sample_id_x, patient_id, sort_short_x, tumor_subsite, therapy) %>% 
  arrange(sample_id_x)

comp_tbl_sample_sort %>% 
  mutate(sort_short_x = cell_sort) %>% 
  unite(sample_id_x, patient_id, sort_short_x, tumor_subsite_rx) %>% 
  select(sample_id_x, cluster_label_sub, n, nrel, nsum) %>% 
  left_join(sample_id_x_tbl, by = "sample_id_x") %>% 
  write_tsv(paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_subtype_compositions.tsv"))

ybreaks <- c(500, 1000, 2000, 4000, 6000, 8000, 10000, 15000, 20000)

max_cells_per_sample <- max(comp_tbl_sample_sort$nsum)
ymaxn <- ybreaks[max_cells_per_sample < ybreaks][1]

comp_plot_wrapper <- function(y = "nrel", switch = NULL) {
  if (y == "nrel") ylab <- paste0("Fraction\nof cells [%]")
  if (y == "n") ylab <- paste0("Number\nof cells")
  p <- ggplot(comp_tbl_sample_sort, 
              aes_string("tumor_subsite_rx", y, fill = "cluster_label_sub")) + 
    facet_grid(~patient_id, space = "free", scales = "free", switch = switch) +
    coord_cartesian(clip = "off") + 
    scale_fill_manual(values = clrs$cluster_label_sub[[coi]]) + 
    theme(axis.text.x = element_blank(),
          axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5, 
                                      margin = margin(0, -0.4, 0, 0, unit = "npc")),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.line.x = element_blank(),
          strip.text.y = element_blank(),
          strip.text.x = element_blank(),
          strip.background.y = element_blank(),
          strip.background.x = element_blank(),
          plot.margin = margin(0, 0, 0, 0, "npc")) + 
    labs(x = "",
         y = ylab) +
    guides(fill = FALSE)
  if (y == "nrel") p <- p + 
    geom_bar(stat = "identity") +
    scale_y_continuous(expand = c(0, 0), 
                       breaks = c(0, 50, 100), 
                       labels = c("0", "50", "100"))
  if (y == "n") p <- p + 
    geom_bar(stat = "identity", position = position_stack()) +
    scale_y_continuous(expand = c(0, 0), 
                       breaks = c(0, ymaxn/2, ymaxn),
                       limits = c(0, ymaxn),
                       labels = c("", ymaxn/2, ymaxn)) +
    expand_limits(y = c(0, ymaxn)) +
    theme(panel.grid.major.y = element_line(linetype = 1, color = "grey90", size = 0.5))
  return(p)
} 

common_label_layers <- list(
  geom_tile(color = "white", size = 0.15),
  theme(axis.text.x = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.line.x = element_blank(),
        strip.text = element_blank(),
        strip.background = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "npc")),
  scale_y_discrete(expand = c(0, 0)),
  labs(y = ""),
  guides(fill = FALSE),
  facet_grid(~patient_id, 
             space = "free", scales = "free")
)

comp_label_site <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_supersite, patient_id), 
       aes(tumor_subsite_rx, label_supersite, 
           fill = tumor_supersite)) + 
  scale_fill_manual(values = clrs$tumor_supersite) +
  common_label_layers

comp_label_rx <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_therapy, consensus_signature, patient_id), 
       aes(tumor_subsite_rx, label_therapy, 
           fill = therapy)) + 
  scale_fill_manual(values = c(`post-Rx` = "gold3", `pre-Rx` = "steelblue")) +
  common_label_layers

comp_label_mutsig <- ggplot(distinct(comp_tbl_sample_sort, tumor_subsite_rx, therapy, tumor_supersite, label_mutsig, consensus_signature, patient_id), 
       aes(tumor_subsite_rx, label_mutsig, 
           fill = consensus_signature)) + 
  scale_fill_manual(values = clrs$consensus_signature) +
  common_label_layers

patient_label_tbl <- distinct(comp_tbl_sample_sort, patient_id, .keep_all = T)

comp_label_patient_id <- ggplot(patient_label_tbl, aes(tumor_subsite_rx, label_patient_id)) + 
  scale_fill_manual(values = clrs$consensus_signature) +
  geom_text(aes(tumor_subsite_rx, label_patient_id, label = patient_id)) +
  facet_grid(~patient_id, 
             space = "free", scales = "free") +
  coord_cartesian(clip = "off") + 
  theme_void() +
  theme(strip.text = element_blank(),
        strip.background = element_blank(),
        plot.margin = margin(0, 0, 0, 0, "npc"))

hist_plot_wrapper <- function(x = "nrel") {
  if (x == "nrel") {
    xlab <- paste0("Fraction of cells [%]")
    bw <- 5
  }
  if (x == "log10n") {
    xlab <- paste0("Number of cells")
    bw <- 0.2
  }
  p <- ggplot(comp_tbl_sample_sort) +
    ggridges::geom_density_ridges(
      aes_string(x, "cluster_label_sub", fill = "cluster_label_sub"), color = "black",
      stat = "binline", binwidth = bw, scale = 3) +
    facet_grid(label_supersite~., 
               space = "free", scales = "free") +
    scale_fill_manual(values = clrs$cluster_label_sub[[coi]]) + 
    theme(axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.title.y = element_blank(),
          axis.line.y = element_blank(),
          strip.text = element_blank(),
          strip.background = element_blank(),
          plot.margin = margin(0, 0, 0, 0, "npc")) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    guides(fill = FALSE) +
    labs(x = xlab)
  if (x == "log10n") p <- p + expand_limits(x = c(0, 3)) + 
    scale_x_continuous(expand = c(0, 0), 
                       labels = function(x) parse(text = paste("10^", x)))
  return(p)
}

pcomp1 <- comp_plot_wrapper("n")
pcomp2 <- comp_plot_wrapper("nrel")

pcomp_grid <- plot_grid(comp_label_patient_id, 
                        pcomp1, pcomp2, 
                        comp_label_site, comp_label_rx, comp_label_mutsig,
                        ncol = 1, align = "v", 
                        rel_heights = c(0.15, 0.33, 0.33, 0.06, 0.06, 0.06))

phist1 <- hist_plot_wrapper("log10n")

pcomp_hist_grid <- ggdraw() +
  draw_plot(pcomp_grid, x = 0.01, y = 0, width = 0.85, height = 1) +
  draw_plot(phist1, x = 0.87, y = 0.05, width = 0.12, height = 0.8)

pcomp_hist_grid

# ggsave(paste0("_fig/02_composition_v6_",coi,".pdf"), pcomp_hist_grid, width = 10, height = 2)

4.3.3 site specific cluster enrichment

comp_tbl_z <- comp_tbl_sample_sort %>% 
  filter(therapy == "pre-Rx",
         !(tumor_supersite %in% c("Ascites", "Other"))) %>% 
  group_by(patient_id, cluster_label_sub) %>% 
  arrange(patient_id, cluster_label_sub, nrel) %>% 
  mutate(rank = row_number(nrel),
         z_rank = scales::rescale(rank)) %>% 
  mutate(mean_nrel = mean(nrel, na.rm = T),
         sd_nrel = sd(nrel, na.rm = T),
         z_nrel = (nrel - mean_nrel) / sd_nrel) %>% 
  ungroup()

ggplot(comp_tbl_z) +
  geom_boxplot(aes(tumor_supersite, z_nrel, color = tumor_supersite), 
               outlier.shape = NA) +
  geom_point(aes(tumor_supersite, z_nrel, color = tumor_supersite), position = position_jitter(), size = 0.1) +
  facet_grid(~cluster_label_sub, scales = "free_x", space = "free_x") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        aspect.ratio = 5) +
  scale_color_manual(values = clrs$tumor_supersite)

ggplot(comp_tbl_z) +
  geom_boxplot(aes(tumor_supersite, z_rank, color = tumor_supersite), 
               outlier.shape = NA) +
  geom_point(aes(tumor_supersite, z_rank, color = tumor_supersite), position = position_jitter(), size = 0.1) +
  facet_grid(~cluster_label_sub, scales = "free_x", space = "free_x") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        aspect.ratio = 5) +
  scale_color_manual(values = clrs$tumor_supersite)

4.4 CD8 T cell diffusion maps

## subset CD8 cells
# seu_obj_cd8 <- seu_obj_sub_sub %>%
#   subset(subset = cluster_label_sub %in% grep("^CD8.T", unique(seu_obj_sub_sub$cluster_label_sub), value = T)) %>%
#   RunUMAP(dims = 1:50, reduction = "harmony", reduction.name = "umapharmony", reduction.key = "umapharmony_")
# 
# write_rds(seu_obj_cd8, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed.rds")

seu_obj_cd8 <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed.rds")

seu_obj_cd8$cluster_label_sub <- setNames(cluster_label_extended[seu_obj_cd8$cluster_extended], seu_obj_cd8$cell_id)

write_rds(seu_obj_cd8, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed.rds")

# root_cell <- "SPECTRUM-OV-041_S1_CD45P_RIGHT_FALLOPIAN_TUBE_ACCCTCATCCAAGCAT"
# seu_obj_cd8_sub <- subset(seu_obj_cd8, cells = c(root_cell, colnames(seu_obj_cd8)[colnames(seu_obj_cd8)!=root_cell][-1]))
# 
# dc_obj <- DiffusionMap(seu_obj_cd8_sub@reductions$harmony@cell.embeddings, k = 100)
# dc_mat <- dc_obj@eigenvectors
# colnames(dc_mat) <- paste0("DC_", 1:ncol(dc_mat))
# seu_obj_cd8_sub[["DC"]] <- CreateDimReducObject(embeddings = dc_mat, key = "DC_", assay = DefaultAssay(seu_obj_cd8_sub))
# 
# write_rds(seu_obj_cd8_sub, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")
# 
# dpt_obj <- DPT(dc_obj)
# write_rds(dpt_obj, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_sub_dpt.rds")
# seu_obj_cd8_sub$DPT1 <- dpt_obj$DPT1
# write_rds(seu_obj_cd8_sub, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")

seu_obj_cd8_sub <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")

seu_obj_cd8_sub$cluster_label_sub <- setNames(cluster_label_extended[seu_obj_cd8_sub$cluster_extended], seu_obj_cd8_sub$cell_id)

write_rds(seu_obj_cd8_sub, "/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/CD8.T_processed_filtered.rds")

# plot_data_sub_sub <- as_tibble(FetchData(seu_obj_sub_sub, c(myfeatures, "cluster_label"))) %>% 
#   left_join(meta_tbl, by = "sample") %>% 
#   mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
#          tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
#   mutate(cell_id = colnames(seu_obj_sub_sub),
#          cluster_label = ordered(cluster_label, levels = my_subtypes),
#          )
#   
# if (cell_sort == "CD45+") {
# plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
# }
# 
# if (cell_sort == "CD45-") {
# plot_data_sub_sub <- filter(plot_data_sub_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
# }
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label)) + 
#   common_layers_disc +
#   ggtitle("Sub cluster") +
#   #facet_wrap(~cluster_label) +
#   scale_color_manual(values = clrs$cluster_label[[coi]])
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
#   common_layers_disc +
#   ggtitle("Patient") +
#   # facet_wrap(~therapy) +
#   scale_color_manual(values = clrs$patient_id_short)
# 
# ggplot(plot_data_sub_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
#   # geom_point(aes(umapharmony_1, umapharmony_2), 
#   #            color = "grey90", size = 0.01, 
#   #            data = select(plot_data_sub_sub, -tumor_supersite)) +
#   common_layers_disc +
#   ggtitle("Site") +
#   # facet_wrap(~therapy) +
#   scale_color_manual(values = clrs$tumor_supersite)
# 
# write_tsv(select(plot_data_sub_sub, cell_id, everything(), -umapharmony_1, -umapharmony_2, -contains("RNA_")), paste0("/work/shah/uhlitzf/data/SPECTRUM/freeze/v6/", coi, "_embedding_sub.tsv"))

4.4.1 UMAP & DC

plot_data_cd8_sub <- as_tibble(FetchData(seu_obj_cd8_sub, c(myfeatures, "cluster_label_sub", "DC_1", "DC_2"))) %>% 
  left_join(meta_tbl, by = "sample") %>% 
  mutate(patient_id = str_remove_all(patient_id, "SPECTRUM-OV-"),
         tumor_supersite = ordered(tumor_supersite, levels = names(clrs$tumor_supersite))) %>% 
  mutate(cell_id = colnames(seu_obj_cd8_sub),
         cluster_label_sub = ordered(cluster_label_sub, levels = names(clrs$cluster_label_sub[[coi]])))
  
if (cell_sort == "CD45+") {
plot_data_cd8_sub <- filter(plot_data_cd8_sub, sort_parameters != "singlet, live, CD45-", !is.na(tumor_supersite))
}

if (cell_sort == "CD45-") {
plot_data_cd8_sub <- filter(plot_data_cd8_sub, sort_parameters != "singlet, live, CD45+", !is.na(tumor_supersite))
}

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label_sub)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label_sub) +
  scale_color_manual(values = clrs$cluster_label_sub[[coi]])

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_cd8_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = cluster_label_sub)) + 
  common_layers_disc +
  ggtitle("Sub cluster") +
  #facet_wrap(~cluster_label_sub) +
  scale_color_manual(values = clrs$cluster_label_sub[[coi]])

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = patient_id_short)) + 
  common_layers_disc +
  ggtitle("Patient") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$patient_id_short)

ggplot(plot_data_cd8_sub, aes(DC_1, DC_2, color = tumor_supersite)) + 
  common_layers_disc +
  ggtitle("Site") +
  # facet_wrap(~therapy) +
  scale_color_manual(values = clrs$tumor_supersite)

5 session info

devtools::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.2 (2019-12-12)
##  os       Debian GNU/Linux 10 (buster)
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Etc/UTC                     
##  date     2021-02-25                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version    date       lib
##  abind                  1.4-5      2016-07-21 [2]
##  ape                    5.3        2019-03-17 [2]
##  assertthat             0.2.1      2019-03-21 [2]
##  backports              1.1.10     2020-09-15 [1]
##  bibtex                 0.4.2.2    2020-01-02 [2]
##  Biobase                2.46.0     2019-10-29 [2]
##  BiocGenerics           0.32.0     2019-10-29 [2]
##  BiocParallel           1.20.1     2019-12-21 [2]
##  bitops                 1.0-6      2013-08-17 [2]
##  boot                   1.3-24     2019-12-20 [3]
##  broom                  0.7.2      2020-10-20 [1]
##  callr                  3.4.2      2020-02-12 [1]
##  car                    3.0-8      2020-05-21 [1]
##  carData                3.0-4      2020-05-22 [1]
##  caTools                1.17.1.4   2020-01-13 [2]
##  cellranger             1.1.0      2016-07-27 [2]
##  class                  7.3-15     2019-01-01 [3]
##  cli                    2.0.2      2020-02-28 [1]
##  cluster                2.1.0      2019-06-19 [3]
##  codetools              0.2-16     2018-12-24 [3]
##  colorblindr          * 0.1.0      2020-01-13 [2]
##  colorspace           * 1.4-2      2019-12-29 [2]
##  cowplot              * 1.0.0      2019-07-11 [2]
##  crayon                 1.3.4      2017-09-16 [1]
##  curl                   4.3        2019-12-02 [2]
##  data.table             1.12.8     2019-12-09 [2]
##  DBI                    1.1.0      2019-12-15 [2]
##  dbplyr                 2.0.0      2020-11-03 [1]
##  DelayedArray           0.12.2     2020-01-06 [2]
##  DEoptimR               1.0-8      2016-11-19 [1]
##  desc                   1.2.0      2018-05-01 [2]
##  destiny              * 3.0.1      2020-01-16 [1]
##  devtools               2.2.1      2019-09-24 [2]
##  digest                 0.6.25     2020-02-23 [1]
##  dplyr                * 1.0.2      2020-08-18 [1]
##  e1071                  1.7-3      2019-11-26 [1]
##  ellipsis               0.3.1      2020-05-15 [1]
##  evaluate               0.14       2019-05-28 [2]
##  fansi                  0.4.1      2020-01-08 [2]
##  farver                 2.0.3      2020-01-16 [1]
##  fitdistrplus           1.0-14     2019-01-23 [2]
##  forcats              * 0.5.0      2020-03-01 [1]
##  foreign                0.8-74     2019-12-26 [3]
##  formattable            0.2.0.1    2016-08-05 [1]
##  fs                     1.5.0      2020-07-31 [1]
##  future                 1.15.1     2019-11-25 [2]
##  future.apply           1.4.0      2020-01-07 [2]
##  gbRd                   0.4-11     2012-10-01 [2]
##  gdata                  2.18.0     2017-06-06 [2]
##  generics               0.0.2      2018-11-29 [2]
##  GenomeInfoDb           1.22.0     2019-10-29 [2]
##  GenomeInfoDbData       1.2.2      2020-01-14 [2]
##  GenomicRanges          1.38.0     2019-10-29 [2]
##  ggplot.multistats      1.0.0      2019-10-28 [1]
##  ggplot2              * 3.3.2      2020-06-19 [1]
##  ggrepel                0.8.1      2019-05-07 [2]
##  ggridges               0.5.2      2020-01-12 [2]
##  ggthemes               4.2.0      2019-05-13 [1]
##  globals                0.12.5     2019-12-07 [2]
##  glue                   1.3.2      2020-03-12 [1]
##  gplots                 3.0.1.2    2020-01-11 [2]
##  gridExtra              2.3        2017-09-09 [2]
##  gtable                 0.3.0      2019-03-25 [2]
##  gtools                 3.8.1      2018-06-26 [2]
##  haven                  2.3.1      2020-06-01 [1]
##  hexbin                 1.28.0     2019-11-11 [2]
##  hms                    0.5.3      2020-01-08 [1]
##  htmltools              0.5.1.1    2021-01-22 [1]
##  htmlwidgets            1.5.1      2019-10-08 [2]
##  httr                   1.4.2      2020-07-20 [1]
##  ica                    1.0-2      2018-05-24 [2]
##  igraph                 1.2.5      2020-03-19 [1]
##  IRanges                2.20.2     2020-01-13 [2]
##  irlba                  2.3.3      2019-02-05 [2]
##  jsonlite               1.7.1      2020-09-07 [1]
##  KernSmooth             2.23-16    2019-10-15 [3]
##  knitr                  1.26       2019-11-12 [2]
##  labeling               0.3        2014-08-23 [2]
##  laeken                 0.5.1      2020-02-05 [1]
##  lattice                0.20-38    2018-11-04 [3]
##  lazyeval               0.2.2      2019-03-15 [2]
##  leiden                 0.3.1      2019-07-23 [2]
##  lifecycle              0.2.0      2020-03-06 [1]
##  listenv                0.8.0      2019-12-05 [2]
##  lmtest                 0.9-37     2019-04-30 [2]
##  lsei                   1.2-0      2017-10-23 [2]
##  lubridate              1.7.9.2    2020-11-13 [1]
##  magrittr             * 2.0.1      2020-11-17 [1]
##  MASS                   7.3-51.5   2019-12-20 [3]
##  Matrix                 1.2-18     2019-11-27 [3]
##  matrixStats            0.56.0     2020-03-13 [1]
##  memoise                1.1.0      2017-04-21 [2]
##  metap                  1.2        2019-12-08 [2]
##  mnormt                 1.5-5      2016-10-15 [2]
##  modelr                 0.1.8      2020-05-19 [1]
##  multcomp               1.4-12     2020-01-10 [2]
##  multtest               2.42.0     2019-10-29 [2]
##  munsell                0.5.0      2018-06-12 [2]
##  mutoss                 0.1-12     2017-12-04 [2]
##  mvtnorm                1.0-12     2020-01-09 [2]
##  nlme                   3.1-143    2019-12-10 [3]
##  nnet                   7.3-12     2016-02-02 [3]
##  npsurv                 0.4-0      2017-10-14 [2]
##  numDeriv               2016.8-1.1 2019-06-06 [2]
##  openxlsx               4.1.5      2020-05-06 [1]
##  pbapply                1.4-2      2019-08-31 [2]
##  pcaMethods             1.78.0     2019-10-29 [2]
##  pillar                 1.4.6      2020-07-10 [1]
##  pkgbuild               1.0.6      2019-10-09 [2]
##  pkgconfig              2.0.3      2019-09-22 [1]
##  pkgload                1.0.2      2018-10-29 [2]
##  plotly                 4.9.1      2019-11-07 [2]
##  plotrix                3.7-7      2019-12-05 [2]
##  plyr                   1.8.5      2019-12-10 [2]
##  png                    0.1-7      2013-12-03 [2]
##  prettyunits            1.1.1      2020-01-24 [1]
##  processx               3.4.2      2020-02-09 [1]
##  progeny              * 1.11.3     2020-10-22 [1]
##  proxy                  0.4-24     2020-04-25 [1]
##  ps                     1.3.2      2020-02-13 [1]
##  purrr                * 0.3.4      2020-04-17 [1]
##  R.methodsS3            1.7.1      2016-02-16 [2]
##  R.oo                   1.23.0     2019-11-03 [2]
##  R.utils                2.9.2      2019-12-08 [2]
##  R6                     2.4.1      2019-11-12 [1]
##  ranger                 0.12.1     2020-01-10 [1]
##  RANN                   2.6.1      2019-01-08 [2]
##  rappdirs               0.3.1      2016-03-28 [2]
##  RColorBrewer           1.1-2      2014-12-07 [2]
##  Rcpp                   1.0.4      2020-03-17 [1]
##  RcppAnnoy              0.0.16     2020-03-08 [1]
##  RcppEigen              0.3.3.7.0  2019-11-16 [2]
##  RcppHNSW               0.2.0      2019-09-20 [2]
##  RcppParallel           4.4.4      2019-09-27 [2]
##  RCurl                  1.98-1.1   2020-01-19 [1]
##  Rdpack                 0.11-1     2019-12-14 [2]
##  readr                * 1.4.0      2020-10-05 [1]
##  readxl               * 1.3.1      2019-03-13 [2]
##  rematch                1.0.1      2016-04-21 [2]
##  remotes                2.1.0      2019-06-24 [2]
##  reprex                 0.3.0      2019-05-16 [2]
##  reshape2               1.4.3      2017-12-11 [2]
##  reticulate             1.14       2019-12-17 [2]
##  rio                    0.5.16     2018-11-26 [1]
##  rlang                  0.4.8      2020-10-08 [1]
##  rmarkdown              2.0        2019-12-12 [2]
##  robustbase             0.93-6     2020-03-23 [1]
##  ROCR                   1.0-7      2015-03-26 [2]
##  rprojroot              1.3-2      2018-01-03 [2]
##  RSpectra               0.16-0     2019-12-01 [2]
##  rstudioapi             0.11       2020-02-07 [1]
##  rsvd                   1.0.3      2020-02-17 [1]
##  Rtsne                  0.15       2018-11-10 [2]
##  rvest                  0.3.6      2020-07-25 [1]
##  S4Vectors              0.24.2     2020-01-13 [2]
##  sandwich               2.5-1      2019-04-06 [2]
##  scales                 1.1.0      2019-11-18 [2]
##  scatterplot3d          0.3-41     2018-03-14 [1]
##  sctransform            0.2.1      2019-12-17 [2]
##  SDMTools               1.1-221.2  2019-11-30 [2]
##  sessioninfo            1.1.1      2018-11-05 [2]
##  Seurat               * 3.1.2      2019-12-12 [2]
##  SingleCellExperiment   1.8.0      2019-10-29 [2]
##  smoother               1.1        2015-04-16 [1]
##  sn                     1.5-4      2019-05-14 [2]
##  sp                     1.4-2      2020-05-20 [1]
##  stringi                1.5.3      2020-09-09 [1]
##  stringr              * 1.4.0      2019-02-10 [1]
##  SummarizedExperiment   1.16.1     2019-12-19 [2]
##  survival               3.1-8      2019-12-03 [3]
##  testthat               2.3.2      2020-03-02 [1]
##  TFisher                0.2.0      2018-03-21 [2]
##  TH.data                1.0-10     2019-01-21 [2]
##  tibble               * 3.0.4      2020-10-12 [1]
##  tidyr                * 1.1.2      2020-08-27 [1]
##  tidyselect             1.1.0      2020-05-11 [1]
##  tidyverse            * 1.3.0      2019-11-21 [2]
##  tsne                   0.1-3      2016-07-15 [2]
##  TTR                    0.23-6     2019-12-15 [1]
##  usethis                1.5.1      2019-07-04 [2]
##  uwot                   0.1.5      2019-12-04 [2]
##  vcd                    1.4-7      2020-04-02 [1]
##  vctrs                  0.3.5      2020-11-17 [1]
##  VIM                    6.0.0      2020-05-08 [1]
##  viridis              * 0.5.1      2018-03-29 [2]
##  viridisLite          * 0.3.0      2018-02-01 [2]
##  withr                  2.3.0      2020-09-22 [1]
##  xfun                   0.12       2020-01-13 [2]
##  xml2                   1.3.2      2020-04-23 [1]
##  xts                    0.12-0     2020-01-19 [1]
##  XVector                0.26.0     2019-10-29 [2]
##  yaml                   2.2.1      2020-02-01 [1]
##  zip                    2.0.4      2019-09-01 [1]
##  zlibbioc               1.32.0     2019-10-29 [2]
##  zoo                    1.8-7      2020-01-10 [2]
##  source                                 
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  Bioconductor                           
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Github (clauswilke/colorblindr@1ac3d4d)
##  R-Forge (R 3.6.2)                      
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  Bioconductor                           
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Github (saezlab/progeny@94bea1c)       
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.3)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.3)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
##  CRAN (R 3.6.2)                         
##  Bioconductor                           
##  CRAN (R 3.6.2)                         
## 
## [1] /home/uhlitzf/R/lib
## [2] /usr/local/lib/R/site-library
## [3] /usr/local/lib/R/library